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Abstract: The ghost free massive gravity modified Friedmann equations at cosmic 
scale and provided an explanation of cosmic acceleration without dark energy. We 
analyzed the cosmological solutions of the massive gravity in detail and confronted 
the cosmological model with current observational data. We found that the model 
parameters and 04 which are the coefficients of the third and fourth order nonlin- 
ear interactions cannot be constrained by current data at the background level. The 
mass of graviton is found to be the order of current Hubble constant if 03 = 04 = 0, 
and the mass of graviton can be as small as possible in the most general case. 
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1 Introduction 

A lot of efforts have been made to understand the accelerating expansion of the 
universe discovered by the observations of Type la supernovae (SNe la) in 1998 [1, 2]. 
Although the economic explanation of the acceleration is a cosmological constant 
which is consistent with all observations, the smallness of the cosmological constant 
and other problems such as the coincidence problem motivated the modification of 
the theory of gravity. One way of modifying gravity is to add a small mass to 
graviton. In fact, Fierz and Pauli made the first attempt to consider a theory of 
gravity with massive graviton [3]. They add a quadratic mass term im?'{h^uh^'^ — 
h?) for linear gravitational perturbations h^i, to the action, which breaks the gauge 
invariance of general relativity. However, the linear theory with the Fierz- Pauli mass 
does not recover general relativity in the massless limit m — >■ 0, which leads to the 
contradiction with solar system tests due to the vDVZ discontinuity [4, 5]. The 
discontinuity can be overcome by introducing nonlinear interactions with the help of 
Vainshtein mechanism [6]. Along this line, Dvali, Gabadadze and Porrati proposed 
a model of massive gravity in the context of extra dimensions which modifies general 
relativity at the cosmological scale and admits a self-accelerating solution with dust 
matter only [7]. 

On the other hand, in the language of Stiickelberg fields, the nonlinear terms 
usually contain more than two time derivatives which present the Bouldware-Deser 
(BD) ghost [8]. To remove the ghost, nonlinear interactions with higher derivatives 
are added order by order in perturbation theory so that they re-sum to be a total 
derivative. Recently, de Rham, Gabadadze and Tolley (dGRT) successfully con- 
structed a nonlinear theory of massive gravity [9] that is free from BD ghost [10]. 
Cosmological solutions with self acceleration for massive gravity were then sought 
by several groups [11-20]. Gumrukcuoglu, Lin, and Mukohyamafound found that a 
de-Sitter solution with an effective cosmological constant proportional to the mass of 
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graviton exists for a spatially open universe in the dGRT model of massive gravity 
[13]. The same solution was then found for spatially fiat universe in [14]. When 
the parameters in the dGRT theory take some particular values, Kobayashi, Siino, 
Yamaguchi and Yoshida found that the solution also existed for a universe with arbi- 
trary spatial curvature [17]. The same solution was obtained by different group with 
different method for some particular case, and they all took the reference metric to 
be Minkowski. Langlois and Naruko took a different approach by assuming the refer- 
ence metric to be de-Sitter and found more general cosmological solutions in addition 
to the cosmological constant solution [20]. These new cosmological solutions opened 
another door to the understanding of cosmic acceleration. Cosmological solutions 
for the ghost-free bi-gravity were also found and confronted with observational data 
[21-25]. 

In this paper, we focus on the cosmological solutions found in [20] for dGRT 
massive gravity [9]. The Friedmann equations are modified so that it is possible to 
explain the cosmic acceleration. We apply the SNLS3 SNe la data [26], the baryon 
acoustic oscillation (BAO) data [27] and the 7-year Wilson Microwave Anisotropy 
Probe (WMAP7) data [28] to constrain the parameters in dGRT massive gravity. 

2 massive gravity 

In this paper, we study the ghost free theory of massive gravity proposed by [9], 

Ml 
2 

where mg is the mass of graviton, the nonlinear higher derivative terms for the 
massive graviton is 

U=U2 + a^Ui + 04^4, (2.2) 
= [IC]' - [IC% (2.3) 
ZY3 = [/C]3-3[/C][/C2]+2[/C3], (2.4) 
= [/C]4 - 6[/C]'[/C'] + 8[/C3][/C] - Q[IC% (2.5) 

and 

/ci; = s^, - {Vex, (2.6) 

The tensor S^^ is defined by four Stiickelberg fields 0° as 

= d^rd,<j)%,. (2.7) 

The reference metric rjab is arbitrary and it is usually taken to be Minkowski. For 
an open universe, Gumrukcuoglu, Lin, and Mukohyamafound found the first cosmo- 
logical solution with an effective cosmological constant proportional to the mass of 
graviton [13], 

Ae// = -ml{X± - 1)[(1 + 3«3)X± - 3(1 + «3)], (2.8) 



8=^-^1 d^xV^iR + mp) + Sm, (2.1) 
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where 

1 + Gas + 12tt4 ± v^l + 3^3 + 9al - 12^4 
- 3(a3 + 4«4) • ^^-^^ 

It is obvious that + Aa^ ^ for this solution, and two branches exist. The same 
solution was then found in [14] for a fiat universe. It is natural to think that this 
solution should exist for a closed universe. If the parameters as and take the 
particular value 

"3=2(0-1), ^4 = — (a^ - a + 1), (2.10) 

the same cosmological constant solution with Ag// = m?g/a independent of the curva- 
ture of the universe was found in [17]. All these results are based on the assumption 
that the reference metric is Minkowski, and the method of obtaining the solution 
cannot be generalized to the other case. 

Langlois and Naruko took a different approach and assumed de Sitter metric for 
the reference metric [20], 

riabd(j)''d(j)^ = -dT^ + hl{T)-iijdX'dX\ (2.11) 

where the Stiickelberg fields are assumed to be 0° = T = fit), = X* = x*, so that 
the tensor S^^, takes the homogeneous and isotropic form, 

S^, = Diag{-/^ &^[/(t)]7^,}, (2.12) 

and the functions hkiT) {k = 0, ±1) are 

6o(T) = e^^^, 6_i(T) = H^^ sinh{H^T), bi{T) = R-^ cosh{H^T). 

Varying the action (2.1) with respect to the lapse function N(t) and scale factor 
a{t), we obtain Friedmann equations 

2H + 3H^ + ^ = -^{p,„+p,), (2.14) 



where the effective energy density pg and pressure Pg for the massive graviton are, 

(2.15) 



Pa = ^J^iMf] - «){6(1 + 2«3 + 2a,)a^ - (3 + 15«3 + 2Aa^)ah[f] 



Pg =^F^{[6 + 12«3 + 12^4 - (3 + 9^3 + 12a4)/]a' 
- 2[3 + 9^3 + 12^4 - (1 + 603 + I2ai)f]ahk[f] 
+ [1 + 6^3 + 12^4 - 3(^3 + 4a4)/]&fc[/]}. 



(2.16) 
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Varying the action (2.1) with respect to the function f{t), we obtain three branches 
of cosmological solutions [20], the first two solutions bk[f(t)] = X± a(t) correspond to 
the effective cosmological constant and are independent of the explicit form of bk{f) 
as long as the function is invertible.^ The third solution is [19, 20] 



df 



a 

N' 



(2.17) 



For the flat case, k = 0, substituting the de Sitter function bo[f{t)] = e^"^^*^ into 
equations (2.17) and (2.15), we obtain the effective energy density and pressure for 
the massive graviton. 



Pg = m^M^i 



-6(1 + 2^3 + 2a4) + 9(1 + 3a3 + 4a4) — 



Jj2 Jj3 

-3(1 + 6^3 + 12a4)-^ + 3(^3 + 4^4)-^ 



(2.18) 



Pg 



-3(1 + 3^3 + 4^4) + 2(1 + 6a3 + 12a4)7r 



J^2 

-3(^3 4^4) — 



(2.19) 



Note that when H{z) = He, Pg = and the contribution to the energy density from 
massive gravity is zero at this moment. For the fiat case, substituting equations 
(2.18) and (2.19) into Friedmann equations (2.13) and (2.14), we get 



Hi H, 



H'^(z) H(z) 
[as + 404)—^ + (1 + 6^3 + 12^4) - 3(1 + 3^3 + 4^4) 



-E^(z)+nm(l+Z 



He 

2 

,3(i+«,„)_2_^(l + 2«3 + 2a4). 



(2.20) 



H, 







H 



TTl 

-2E\z) + -^E(z) 



3(1 + 3^3 + 4^4)77^ - 2(1 + 6^3 + l2ai)E{z) 
Ho 



^0 



+'i{a^ + Aa^)-^E\z) 



(2.21) 



3f]„,(l + w;J(l + z)3(i+'"'-), 



where E{z) = H{z)/Hq. The effective equation of state Wg = Pg/pg for the massive 
graviton is 



w„ 



2E\z)H/H^ + Samjl + wjjl + ^)3(i+^^) 

3[E2(z) -fi^(l + Z)3(l+«'m)] 



(2.22) 



Without loss of generality, we assume that 



-/3li/o^ and H^ = P2H0. Note 



that the mass appears in the action as a potential term, so the sign of m can be 



^The effective cosmological constant (2.8) is obtained by substituting the solution hk[f{t)] 
X± a{t) into the energy density of massive graviton (2.15) 
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absorbed into the sign convention of the potential. At the present time z 
E(0) = 1, equation (2.20) gives 



(1 - 



^ 2(1 + 2^3 + 2a4)/3| - 3(1 + Sag + 4a4)/3| + (1 + 603 + \1ai)^2 - («3 + 4^4) ' 

(2.23) 

If 03 + 404 = 0, the cubic term of Hubble parameter in equation (2.20) is absent and 
the cosmological evolution becomes simpler. Therefore we consider this special case 
first. For the special case 4a4 = —03, Friedmann equation (2.20) becomes 



1- §(1 + 803) 
P2 



E\z) + 3(1 + 2a^)^E{z) = n^l + z)3(i+-'") + /3,(2 + 803), 

P2 



with 



/3i 



(2.24) 



(2.25) 



[(2 + 3a3)/32-l-3a3](/32-l) 
When > 1, /3i ^ (1 - f^m)/(2 + 803) and the standard ACDM model with 
cosmological constant Q\ = 1 — is recovered. Note that when ^ 1, the model 
just weakly depends on the parameter as and ^ when 0:3 ^ 1. The deceleration 
parameter in this model is 

q{z) = -1 + 



2-2(l + 3a3)| , , 



3il + 2as)^E-\zy 



(2.26) 



In this case, we have three model parameters ^2^, (^2 and a^. Apparently, the coeffi- 
cient of H'^{z) should be positive, so the model parameters must satisfy the following 
condition 

[(2 + 3a3)/32-l-3a3](/32-l) ' ^ ' ' 

At early times, E{z) ^ 1, the square term E'^{z) dominates the left hand side of 
equation (2.24), the standard cosmology with an effective cosmological constant is 
recovered and the effective matter density is ^^/[l — /3i(l + 3a;3)//3|] instead of fi^. 
To guarantee that equation (2.24) always has solutions, we require that 



A = 9(l + 2«3)'^ + 4 



PI 



/3 



Ul + 3a3 



C{z) > 0, 



(2.28) 



where C{z) = ^^(1 + 2)3(1+"'^) _^ ^^(2 + 3^3). Explicitly, the dimensionless Hubble 
parameter E{z) is 

-3(l + 2a3)/3i//32 + 



E(z) 



2[1 - (1 + 3as)/3,/(3l] 



(2.29) 
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Finally, to enure that E{0) = 1, we require 



- 3(1 + 2as)^ < 2 - 2(1 + 3«3)4. (2.30) 

To better understand the dynamics, we analyze the simplest case 03 = 04 = in 
more details, in which we have only two model parameters and /32- The Hubble 
parameter evolves as 



= ^ ^ on \ /.2^ • (2-3i; 



The condition (2.27) is reduced to 

(1-a 



(2/32-l)(/32 



< 1. (2.32) 



(1) If 1/2 < j32 < 1, then /3i < 0, > and the above condition (2.32) is 
automatically satisfied. Combining this result with condition (2.28), we get 1/2 < 
/32 < (3 - yrU) /4 or (3 + y/H^) /4 < /^a < 1 and > 2H^. Since equation (2.24) is 
a quadratic equation, there are two solutions, and we need to take the solution which 
E[z) increases as the redshift z increases and we require that E{z = 0) = 1. To satisfy 
these requirements, the parameter f32 must be in the region 1/2 < (32 < (3 — A/^m)/4. 

(2) If (32 > 1, then (3i > and < 0. The conditions (2.32) and (2.28) require that 
/32 > (3 + V9-8aj/4. When (32 > 1, A ^ (1 - nj/2, ^ -(1 - nm)H^/2 
and the standard ACDM model is recovered. (3) li < (32 < 1/2, the conditions 
(2.32) and (2.28) require that < /^a < (3 - V9 - 8fi^)/4. When P2 1, (3i ^ 
(1 — Qrn)P2 1 and ^ 0. However, the standard ACDM model can not be 
recovered at early times. Therefore, we don't consider this case. 

Now we consider the general case with + 4a4 7^ 0. Since Pg = when 
H[z) = He, so if He = Hq, we find that VLm = 1 which is inconsistent with current 
observations, so Hc^ Hq. If He < Hq, then in the past z ^ 1, the cubic term H^{z) 
dominates over the quadratic H'^{z) and the linear H{z) terms, so we cannot recover 
the standard cosmology H^ ~ p unless we fine tune the value of m^/i^Q to be very 
small. Therefore, we require (32 = H^/Hq > 1. From equation (2.20), we see that 
the standard cosmology is recovered when Hq < H{z) < H^ ii (32 ^ 1- At very early 
times when H{z) ^ He, the universe evolves according to H^ ~ p. During radiation 
dominated era, the universe evolves faster according to a{t) ~ t^^^ instead of t^/^. 

If the parameters 03 and take the particular values in equation (2.10), then 
Friedmann equation (2.20) becomes 



3^,^1(1 + afE{z) + 'i(32[(3l - /3ia(l + a)]E''{z) + (3ia''E%z) 



(2.33) 
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with 

« = 3(1 - .2 341 

(a2 + 3a + 3)/3| - 3{a + If (31 + 3a(a + 1)132 - a^' 

The deceleration parameter is 

' Pi[aE{z)-{l + a)[32?E{z) + 2l3lE^{zy ' 

When ^ 1, for most values of a, ~ (1 — is negligible, therefore 

the cubic term E^{z) can be neglected, which means that the model is not sensitive 
to the parameter a. The same is true for the most general case, so we expect that 
the model parameters 03 and 0:4 are not well constrained by the observational data 
at the background level. 



3 Observational constraints 

To find out the parameters which are consistent with observatioal dat, we use the 
SNLS3 SNe la [26], BAO [27] and WMAP7 data [28]. The SNLS3 SNe la data 
consists of 123 low-redshift SNe la data with 2; < 0.1 mainly from Calan/Tololo, CfAI, 
CfAII, CfAIII and CSP, 242 SNe la over the redshift range 0.08 <z < 1.06 observed 
from the SNLS [26], 93 intermediate-redshift SNe la data with 0.06 ^ z < 0.4 
observed during the first season of Sloan Digital Sky Survey (SDSS)-II supernova 
(SN) survey [29], and 14 high-redshift SNe la data with z > 0.8 from Hubble Space 
Telescope [30]. To fit the SNL3 data, we need two additional nuisance parameters 
a and /3 in addition to the model parameters p. For the BAO data, we use the 
measurements of the distance ratio dz = rs{zd) / Dv{z) at the redshift z = 0.106 from 
the 6dFGS [31], the measurements of at two redshifts z = 0.2 and z = 0.35 from 
the distribution of galaxies [32], and the measurements of the acoustic parameter 
A{z) = Dv{z)^VtmHl/z at three redshifts z = 0.44, z = 0.6 and z = 0.73 from 
the WiggleZ dark energy survey [27]. In the fitting of BAO data, we need to use 
two more nuisance parameters Qhh"^ and flmh"^- For the WMAP7 data [28], we use 
the measurements of the three derived quantities: the shift parameter R{z*), the 
acoustic index Ia{z*) and the recombination redshift z*. The nuisance parameters 
Qbh"^ and Qmh"^ are again needed when we employ the WMAP7 data. We apply the 
Monte Carlo Markov Chain code [33, 34] to find out the best fitting parameters. 

As we discussed above, for the case + 4^4 = 0, when ^ 1, the model is 
almost equivalent to ACDM model and the model is not sensitive to the values of 
and (32, so the value of (^2 is not bounded from above. By fitting the model to the 
observational data, we find that the minimum value of for /32 < 1 is much bigger 
than that for ^2 > 1. Therefore, the parameter space with /32 < 1 can be ignored 
and we choose (32 > I. Because (32 is unbounded from above, we cut (32 at e^°. 
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When 03 = 04 = 0, fitting the model to the observational data, we get the 
minimum value of = 421.65, the best fit values Qm = 0.27 and /32 = 2.44 and 
P2 > 1-93 at 3a. Note that for the curved ACDM model, the minimum value of 
is 423.98. The marginalized probability distributions of the model parameters 
Qm, In /32 and Pi along with the marginalized contours of /3i and ln/32 are shown 
in Fig. 1. As expected, we get a long flat tail for large (32 and the most probable 
marginalized value of /3i is around the asymptotical value Pi ~ 0.37, so the mass of 
graviton is around 0.6Hq. From the mean likelihood distribution (the dotted line), 
we see that Pi is peaked at its best fitted value. Because P2 is not Gaussian, so the 
derived quantity Pi is not Gaussian either. On the other hand, P2 is highly peaked 
at its best fit value. To see this point clearly, we plot the function of versus P2 
in Fig. 2 by fixing the other parameters at their best fit values. From the plot in 
Fig. 2, we see that the probability of P2 is negligible when it is away from its best 
fit value. Note that this is different from the marginalized probability distribution 
because we neglect the degeneracies between P2 and other parameters including the 
nuisance parameters Qbh'^ and Hq. Due to the two effects discussed above, the peak 
value of Pi in the marginalized likelihood distribution is different from that in the 
mean likelihood distribution. Using the best fit values, we plot the evolutions of the 
deceleration parameter q{z) and the matter density ^^(z) = 87iGpm/{3H'^) in Fig. 
2. 

For the case a^ + Aa^ = 0, we find that the best fit values are Qm = 0.27, P2 = 3.1 
and as = 0.95 with = 421.57. The marginalized probability distributions of the 
model parameters flrn, ln/32, 03 and Pi are shown in Fig. 3. As expected, the model 
weakly depends on the value of 0:3, and nonnegative values are slightly more probable 
than negative values. From the marginalized probability distribution of Pi, we see 
that the most probable value is /3i ~ which means the graviton is almost massless. 
This can be understood from equation (2.25), when P2 ^ 1 and jasl ^ 1, /3i ~ 0. 
For the same reason of non-Gaussian distribution, in the mean likelihood distribution 
(the dotted line). Pi is peaked at its best fitted value. With the best fit values, we 
construct the evolutions of the deceleration parameter q{z) and the matter density 
^lm{z) and they are shown in Fig. 2. 

As discussed above, the parameters 03 and 0:4 are uncorrelated with Qm and 
P2 for the general case, so we donot expect that the parameters 03 and 04 can be 
constrained by the observational data, and the results are similar to that of the 
special case 03 = 04 = 0. If the parameters and take the particular value 
(2.10), the best fitting values are fi^ = 0.28, ln/32 = 14.42 and a = -71.3 with 

= 424.28. Since /32 ^ 1, the probability distributions of P2 and a are almost 
fiat as expected. For the most general case, the best fitting values are Q^n = 0.28, 
ln/32 = 15.83, as = -475.4 and = -362.1 with = 424.28. The fitting resuh 
for general 03 and 04 is almost the same as that when 0^3 and 0:4 take the particular 
values, and as and a^ show a fiat distribution, the plots are shown. Taking into 
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Figure 1. The marginalized probability distributions of the model parameters i^m, ^^(^2 
and /?! and the marginalized contours of /3i and ln/32. The dotted line in the lower left 
panel shows the mean likelihood distribution. 

account the results obtained for the special case, we see that the special case with 
as + 4^4 = fits the results better. 

4 Conclusions 

The cosmological constant solution (2.8) for different case was found by different 
group with different method by assuming the reference metric to be Minkowski. 
Therefore the solution (2.8) should be true in the most general case with arbitrary 
spatial curvature and any model parameters. In fact, the solution was found to be 
true in the most general case for an isotropic and homogeneous universe by taking 
the ansatz (2.12) for the tensor S^j,. Furthermore, new cosmological solution which 
modified Frriemann equation was also found in this approach. Therefore, the new 
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Figure 2. Left panel: versus ln/32 for the special case 03 = 0:4. The other parameters 
are fixed at their best fit values. Right panel: The evolutions of q{z) (black lines) and 
^miz) (blue lines), the solid lines are for the model with 03 = 04 = and the broken lines 
are the model with 0:3 + 4q;4 = 0. 

solution should also be a general solution and more cosmological solutions can be 
found following this approach [12]. 

Fitting the model to the observational data, we find that the best fit values 
^Irn = 0.27 and (32 = 2.4 with = 421.65 when = 04 = 0. We also find that 
(32 > 1.93 at 3cr level and the mass of graviton is around O.GHq. For the special case 
a-s + 4^4 = 0, we find that the best fit values flm = 0.27, (32 = 3.1 and = 0.95 
with = 421.57, which is almost the same as the case with ^3 = 0. In fact, we find 
that as is almost uncorrelated with the parameters Qm and ^2- For the most general 
case, 03 and are uncorrelated with the parameters Qm and P2 and therefore are 
not constrained by current observational data at least at the background level. The 
simple case with + 4^4 = is slightly favored by the observational data, so for 
phenomeno logical interest, we can consider the simple case only. 

Although the model fits the observation as well as the standard ACDM model 
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Figure 3. The marginalized probability distributions of the model parameters Vtm, ln/32, 
a-i and The dotted line in the lower right panel shows the mean likelihood distribution. 



does, the phenomenology of the model is distinctively different from dark energy 
models. As seen in Fig. 2, the matter density exceeds the critical density around 
redshift z = 2 which may make the model testable, and the period of this over-density 
is longer for the model with more degrees of freedom [a^ ^ 0). The consequence of 
the feature and how to detect massive gravity from astrophysical observations need 
to be further studied. 
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